1. Wstęp
  2. Opis metod prognozowania
    1. Modele regresyjne
    2. Prophet
    3. TBATS
    4. Ostateczna predykcja
  3. Przewidywanie
  4. Podsumowanie

1. Wstęp

Wyzwanie konkursowe polega na predykcji wskaźnika bezrobocia na lata 2022-2023, na podstawie danych z 2000-2021. Analizę będziemy przeprowadzać przy użyciu języka R w programie RStudio, do obliczeń i przewidywań użyjemy m.in. z bibliotek takich jak stats, tseries, forecast, prophet, do wizualizacji - ggplot2, a do utworzenia ostatecznego raportu posłużymy się oprogramowaniem RMarkdown.

Wejściowe dane przyjmują formę szeregu czasowego - realizacji procesu stochastycznego, którego dziedziną jest czas - w tym przypadku po 12 miesięcy z 22 lat. Procesem stochastycznym \((X_t)_{t \in T}\) nazywamy rodzinę zmiennych losowych z pewnej przestrzeni probabilistycznej, przyjmującą wartości z przestrzeni mierzalnej.

Dane wskaźnika bezrobocia w latach 2000-2021 przedstawiają się w następujący sposób:

Z podziałem na lata:

Od razu zauważamy, że dane podlegają dostrzegalnemu, nieliniowemu trendowi - wartość wskaźnika spada na początku roku, w okolicy połowy roku wzrasta, by potem delikatnie spadać, lub utrzymywać się aż do października, a na koniec roku znowu wzrasta. Występuje oczywiście sezonowość o okresie 12 miesięcy

By sprawdzić podejrzenia wynikające z wizualnej obserwacji powyższych wykresów, przeprowadzamy dekompozycję STL (Seasonal Decomposition of Time Series by Loess). Ma ona wskazać składniki szeregu czasowego - jego trend, sezonowość oraz reszty.

library(forecast)

ts_stl <- ts(df$Wskaznik, frequency = 12, start = c(2000,12))
autoplot(mstl(ts_stl),col=ts_stl)

Powyższa metoda opiera się na przedstawieniu punktów szeregu czasowego (\(y_i, i \in T\)) jako suma komponentów sezonowości \(s_i\), trendu \(t_i\) i reszty \(r_i\): \[y_i = s_i + t_i + r_i\] oraz estymacja owych komponentów. [1]

Obserwując surowe dane, widzimy pewną anomalię w roku 2020 - wyraźny skok wartości wskaźnika bezrobocia w kwietniu w wyniku wybuchu pandemii COVID-19. To zaburzenie w danych w większości przypadków obniży jakość prognozy, ponieważ trend w tym okresie zostaje naruszony. Z tym problemem możemy poradzić sobie na kilka sposobów. Istnieje opcja podstawienia średniej wartości wskaźnika z każdego miesiąca do odpowiednich miesięcy z 2020. Można wziąć średnie globalne, lub jedynie z kilku ostatnich lat. Nie ma również większych przeszkód, by zupełnie pominąć ten rok w obliczeniach. Zbadamy także ideę, by użyć danych z lat 2000-2019, by “przewidzieć” wartości z 2020, a następnie dokonywać obliczeń przy użyciu nowych danych z 2020 do predykcji 2022-2023. Dokonamy analizy tych metod w rozdziale 3.

2. Opis metod prognozowania

2.1 Modele regresyjne

Jako pierwszą metodę predykcji wybraliśmy regresję liniową ze względu na miesiące. Wartości wskaźnika z n-tego miesiąca z lat 2000-2021 wyznaczają przewidywaną wartość wskaźnika z n-tego miesiąca na lata 2022 i 2023. Ze względu na podatność regresji liniowej na obserwacje odstające zdecydowaliśmy się usunąć rok 2020, ponieważ metoda ta nie wymaga ciągłości danych. [CITATION2]

library(stats)
library(tseries)
library(forecast)

predict_monthly <- function(data, month) {
  monthly_data <- data %>% filter(M.c == month)
  model <- lm(Wskaznik ~ Rok, data = monthly_data)
  future_years <- data.frame(Rok = c(2022, 2023), M.c = month)
  predictionsWith2020 <- predict(model, newdata = future_years)
  return(data.frame(Rok = future_years$Rok, M.c = future_years$M.c, Wskaznik = predictionsWith2020))
}

predictionsWithout2020 <- lapply(1:12, function(m) predict_monthly(data, m))

predictionsWithout2020_df <- do.call(rbind, predictionsWithout2020)
combined_data <- bind_rows(data, predictionsWithout2020_df)

Ponieważ regresja liniowa jest ogólnym modelem, który nie jest dedykowany dla szeregów czasowych, zdecydowaliśmy się rozszerzyć te metodę o model ARIMA (Auto Regressive Integrated Moving Average). [CITATION3]

W modelu tym, część autoregresyjna (AR) jest zasadniczo formą regresji liniowej, gdzie bieżąca wartość szeregu czasowego jest modelowana jako liniowa kombinacja jego przeszłych wartości.

Ponieważ ARIMA również jest podatna na obserwacje odstające oraz wymaga ciągłości danych, zdecydowaliśmy się zastąpić rok 2020 predykcją modelu ARIMA na podstawie lat 2000-2019 i wykorzystać te dane aby przewidzieć wartości wskaźnika na podstawie lat 2000-2021 ze sztucznymi danymi z roku 2020. Mając swiadomość możliwości wystąpienia efektu kaskadowania błędów, porównamy to rozwiązanie z klasycznym zastąpieniem roku 2020 średnimi z pozostałych lat.

Tak wyglądają wartości z 2020, których użyjemy w predykcji lat 2022-2023:

fit_to_2019 <- auto.arima(ts_data)
forecast_2020 <- forecast(fit_to_2019, h=12)

Korzystając z funkcji adf.test sprawdzamy, czy spełniona jest stacjonarność modelowanych danych, która jest kluczowym założeniem modelu ARIMA.

adf.test(ts_data_doubleARIMA)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_data_doubleARIMA
## Dickey-Fuller = -9.5335, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
adf.test(ts_data_mean)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_data_mean
## Dickey-Fuller = -9.4633, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

KOMENTARZ JAKIE WYNIKI TESTOW

Zbadajmy następnie błędy ?????? co na skladowych wykresu

fit_doubleARIMA <- auto.arima(ts_data_doubleARIMA)
forecast_values_doubleARIMA <- forecast(fit_doubleARIMA, h=24)
fit_mean <- auto.arima(ts_data_mean)
forecast_values_mean <- forecast(fit_mean, h=24)

checkresiduals(fit_mean)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1)(0,1,2)[12]
## Q* = 27.002, df = 20, p-value = 0.1352
## 
## Model df: 4.   Total lags used: 24
checkresiduals(fit_doubleARIMA)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1)(0,1,2)[12]
## Q* = 27.762, df = 20, p-value = 0.1152
## 
## Model df: 4.   Total lags used: 24

Na powyższych wykresach błędów modelu widać, że zastosowanie podwójnej predykcji (drugi przypadek) jest teoretycznie lepsze niż model biorący rok 2020 jako średnią z pozostałych lat DLACZEGO. Jednakze, błędy w 2020 są praktycznie zerowe - wynika to z faktu, że rok ten jest dopasowany przez ten sam model. Mimo to zdecydowaliśmy się przyjąć wyniki z modelu korzystającego z podwójnej predykcji. LJUNG BOX TEST??

forecast_values_doubleARIMA <- data.frame(
  Date = seq(as.Date("2022-01-01"), by = "month", length.out = 24),
  Point_Forecast = forecast_values_doubleARIMA$mean[1:24]
)

forecast_values_mean <- data.frame(
  Date = seq(as.Date("2022-01-01"), by = "month", length.out = 24),
  Point_Forecast = forecast_values_mean$mean[1:24]
)

df_plot$Type <- "Actual"
forecast_values_doubleARIMA$Type <- "Double ARIMA"
names(forecast_values_doubleARIMA)[2] <- "Wskaznik"

forecast_values_mean$Type <- "Mean"
names(forecast_values_mean)[2] <- "Wskaznik"

combined_df <- rbind(df_plot[, c("Date", "Wskaznik", "Type")], 
                     forecast_values_doubleARIMA, 
                     forecast_values_mean)

Ostatecznie, dane przewidziane przy użyciu omówionych metod modeli regresyjnych przedstawiają się następująco:

ARIMA <- combined_df[combined_df$Type == 'Double ARIMA',]

2.2 Prophet

Kolejną metodą przewidywania szeregu czasowego jest Prophet, przedstawiony przez Facebooka.[4] Na wyjściowym szeregu dokonujemy dekompozycji w następujący sposób: \[y(t) = g(t) + s(t) + h(t) + e_t\] W tym przypadku \(g(t)\) jest funkcją trendu reprezentującą nieokresowe zmiany wartości szeregu czasowego, \(s(t)\) jest funkcją zmian okresowych (np. miesięcznych), a \(h(t)\) to efekt świąt, występujących w nieregularnych odstępach czasu. Zakładamy, że błąd \(e_t\) ma rozkład normalny. funkcje \(g(t)\), \(s(t)\), \(h(t)\) są estymowane przy użyciu m. in. logistycznego modelu wzrostu oraz szeregów Fouriera, ściślej opisane w [4].

Podstawiając nasze dane uzyskujemy następujące wyniki predykcji:

library(prophet)

df_prophet <- df
df_prophet$date <- as.Date(paste(df$Rok, df$M.c, "01", sep = "-"))
df_prophet <- subset(df_prophet, select = c(date, Wskaznik))
colnames(df_prophet) <- c('ds', 'y')
model_prophet <- prophet(df_prophet)
future <- make_future_dataframe(model_prophet,
                                periods = 24,
                                freq = 'month')
forecast <- predict(model_prophet, future)
forecast[c('ds', 'yhat')] %>%
  tail(n = 24)
##             ds      yhat
## 265 2022-01-01 105.59152
## 266 2022-02-01 100.91219
## 267 2022-03-01  98.81112
## 268 2022-04-01  97.07026
## 269 2022-05-01  97.14200
## 270 2022-06-01  97.91921
## 271 2022-07-01  99.66101
## 272 2022-08-01  99.95366
## 273 2022-09-01  99.73869
## 274 2022-10-01  99.50115
## 275 2022-11-01 101.79561
## 276 2022-12-01 103.00231
## 277 2023-01-01 105.79901
## 278 2023-02-01 101.15992
## 279 2023-03-01  99.35893
## 280 2023-04-01  97.09987
## 281 2023-05-01  97.20263
## 282 2023-06-01  98.16253
## 283 2023-07-01 100.04565
## 284 2023-08-01 100.28066
## 285 2023-09-01  99.94061
## 286 2023-10-01  99.66011
## 287 2023-11-01 101.96982
## 288 2023-12-01 103.25350
prophet_plot_components(model_prophet, forecast)

plot(model_prophet, forecast)

Na pierwszym wykresie widzimy krzywą trendu wyznaczoną przez model, łącznie z przewidzianymi ostatnimi latami wraz z przedziałem ufności (niebieskie pole na końcu krzywej), na drugim rysunku - uśrednione, ogólne roczne zmiany wskaźnika. Ostatni wykres przedstawia dodatkowo porównanie rzeczywistych danych (czarne punkty) z tymi estymowanymi przez Prophet (niebieska linia, wraz z przedziałem ufności).

Zbadajmy teraz jakość predykcji, kiedy pominiemy rok 2020 w rozważaniach:

forecast2[c('ds', 'yhat')] %>%
  tail(n = 24)
##             ds      yhat
## 253 2022-01-01 104.48811
## 254 2022-02-01  99.89273
## 255 2022-03-01  97.55734
## 256 2022-04-01  95.63517
## 257 2022-05-01  95.70601
## 258 2022-06-01  96.53879
## 259 2022-07-01  98.34170
## 260 2022-08-01  98.68818
## 261 2022-09-01  98.51341
## 262 2022-10-01  98.31010
## 263 2022-11-01 100.63690
## 264 2022-12-01 101.82312
## 265 2023-01-01 104.36696
## 266 2023-02-01  99.81617
## 267 2023-03-01  97.72722
## 268 2023-04-01  95.69684
## 269 2023-05-01  95.67828
## 270 2023-06-01  96.49088
## 271 2023-07-01  98.35158
## 272 2023-08-01  98.70011
## 273 2023-09-01  98.48975
## 274 2023-10-01  98.33266
## 275 2023-11-01 100.65959
## 276 2023-12-01 101.92684
prophet_plot_components(model_prophet2, forecast2)

Obserwujemy, że po wyrzuceniu 2020, linia trendu zdecydowanie się wypłaszcza na koniec badanego okresu. Przez to, że model bierze pod uwagę współczynnik świąt, czyli nieregularnych skoków badanej zmiennej, w tym przypadku warto zostawić dane z 2020, zatem ostatecznie bierzemy pierwotną wersję prognozy.

PROPHET <- forecast['yhat'] %>%
  tail(n = 24)

2.3 TBATS

Metoda TBATS należy do popularnej grupy modeli statystycznych - modeli wygładzania wykładniczego. Do tej samej rodziny należy także STL, którą używaliścy do dekompozycji naszych danych.

Najczęściej używany model sezonowy przedstawia się w następujący sposób: [5]

\[y_t = l_{t-1} + b_{t-1} + s_t + d_t\] \[l_t = l_{t-1} + b_{t-1} + \alpha d_t\] \[b_t = b_{t-1} + \beta d_t\] \[s_t = s_{t-m} + \gamma d_t\] gdzie \(m\) to okres cykli sezonowych, \(d_t\) reprezentuje losowy szum w danych, \(l_t\), \(b_t\), \(s_t\) to komponenty odpowiednio: poziomu, trendu i sezonowości szeregu czasowego. Wartości \(\alpha\), \(\beta\) i \(\gamma\) są tak zwanymi parametrami wygładzającymi, a \(l_0\), \(b_0\), \(\{s_{1-m}, ..., s_0\}\) to zmienne wyjściowe. W tym modelu estymuje się właśnie zmienne wyjściowe oraz parametry. Warto dodać, że można estymować 2 składniki sezonowości \(s_t^{(1)}\), \(s_t^{(2)}\) wraz z parametrami \(\gamma_1\) i \(\gamma_2\), jednak w naszym przypadku bierzemy tylko jeden składnik - miesięczny.

Następnie dokonujemy modyfikacji tego modelu, przy uwzględnieniu transformacji Box-Cox, błędów modelu ARMA oraz T wzorców sezonowości. Na koniec, składniki sezonowości przedstawiamy jako ich trygonometryczną reprezentację opartą na szeregach Fouriera. Ze wszystkich składowych modelu powstała jego nazwa: T - trigonometrical, B - Box-Cox transformation, A - ARMA errors, TS - T seasonal patterns.

W celu ewaluacji jakości prognozy, będziemy porównywać ostatnie dane 2 lata z przewidzianymi wartościami na następne 2 lata, korzystając z dwóch metryk:

  1. średni bezwzględny błąd procentowy - MAPE (Mean Absolute Percentage Error): \[\frac{1}{n} \sum_{t=1}^n |\frac{Y_t - P_t}{Y_t}| * 100\%\]
  2. pierwiastek błędu średniokwadratowego - RMSE (Root Mean Squared Error): \[\sqrt{\frac{\sum_{t=1}^n (Y_t - P_t)^2}{n}}\] gdzie \(Y_t\) - rzeczywista wartość, \(P_t\) - prognozowana wartość.

Po podstawieniu naszych danych otrzymujemy następujące wyniki:

ts_tbats <- msts(df$Wskaznik, seasonal.periods = 12)

model.ts_tbats <- tbats(ts_tbats)

model.ts_tbats.forecast <- forecast(model.ts_tbats, h = 24)

plot(model.ts_tbats.forecast, main = 'Prognoza TBATS',
     ylab = 'Wskaznik')

model.ts_tbats.forecast$mean
##          Jan       Feb       Mar       Apr       May       Jun       Jul
## 23 102.29247  98.84214  96.22711  95.65658  95.33452  96.74659  97.75088
## 24 104.10502 100.24764  97.32554  96.53332  96.03627  97.31860  98.21517
##          Aug       Sep       Oct       Nov       Dec
## 23  98.91372  98.21906  98.66880 100.31739 102.23146
## 24  99.29118  98.52024  98.91193 100.51604 102.39415

W wykresie prognozy, niebieska linia oznacza przewidziane dane, wraz z przedziałami ufności na poziomie istotności \(0.2\) oraz \(0.05\). Przewidziane wartości zostały wypisane w tabeli pod nim.

df.test <- tail(df$Wskaznik, n = 12 * 2)

ts_tbats.predict <- predict(model.ts_tbats.forecast, df.test)

ts_tbats.test_vs_predicted <- data.frame(df.test, model.ts_tbats.forecast$mean)
matplot(ts_tbats.test_vs_predicted, type = 'l', lty = 1:2, col = 1:2, ylab = 'Przypadek nr 1')

Na tym rysunku widzimy porównanie wartości z ostatnich dwóch lat (czarna linia) z przewidzianymi wartościami na następne 2 lata (czerwona przerywana linia). Obserwujemy tutaj zauważalne niedopasowanie krzywych na początku okresu - jest to następstwo wzięcia pod uwagę odstających wartości z roku 2020.

MAPE.error(ts_tbats.test_vs_predicted$df.test,
           ts_tbats.test_vs_predicted$model.ts_tbats.forecast.mean)
## [1] 2.085391
RMSE.error(ts_tbats.test_vs_predicted$df.test,
           ts_tbats.test_vs_predicted$model.ts_tbats.forecast.mean)
## [1] 6.855342

Ostatecznie, liczymy wartości błędów MAPE i RMSE - są one zawyżone z wyżej wymienionego powodu.

Aby zmiejszyć wartości błędów, dokonujemy analizy bez 2020:

## MAPE:  1.168697
## RMSE:  3.977294

W przypadku wyrzucenia danych z 2020, zgodnie z oczekiwaniami, widzimy znaczny spadek błędów MAPE i RMSE. Zwróćmy uwagę, że do ewaluacji prognozy wzięliśmy lata 2019 i 2021 zamiast 2020-2021. Sprawdzimy jeszcze predykcję dla przypadków, gdy zamiast wartości z tego roku podstawimy średnią globalną wartość wskaźnika:

## MAPE:  1.009071
## RMSE:  2.643147

Po raz kolejny, wartości błędu zmiejszyły się. Na koniec zbadajmy zachowanie błędów, gdy zamiast 2020 weźmiemy średni wskaźnik z ostatnich 5 lat:

## MAPE:  1.06453
## RMSE:  3.255959

Najmniejsze wartości błędów otrzymujemy w przypadku podstawienia globalnej średniej za wartości z odstającego roku. Decydujemy jednak uwzględnić ostatni przypadek w ostatecznej predykcji, gdyż nie chcemy brać takich wyników, które będą zbyt zbliżone do wartości z ostatnich dwóch lat, ponieważ po raz kolejny, trend byłby zbyt spłaszczony na końcu okresu. Ostatecznie:

TBATS <- model.ts_tbats4.forecast$mean %>% as.data.frame()

2.4 Ostateczna predykcja

Wszystkie wyniki wygenerowane przez omówione modele predykcji przedstawiają się w następujący sposób:

##          Date     ARIMA   Prophet     TBATS
## 1  2022-01-01 102.82608 105.59152 102.15098
## 2  2022-02-01  98.62021 100.91219  98.84227
## 3  2022-03-01  96.13074  98.81112  96.15523
## 4  2022-04-01  95.42150  97.07026  95.20574
## 5  2022-05-01  95.86357  97.14200  94.80567
## 6  2022-06-01  95.77652  97.91921  96.32381
## 7  2022-07-01  97.74766  99.66101  97.41151
## 8  2022-08-01  98.41741  99.95366  98.46648
## 9  2022-09-01  97.32115  99.73869  97.71175
## 10 2022-10-01  97.46430  99.50115  98.19079
## 11 2022-11-01  99.36107 101.79561  99.75118
## 12 2022-12-01 100.32231 103.00231 101.52969
## 13 2023-01-01 103.52965 105.79901 103.32886
## 14 2023-02-01  99.01594 101.15992  99.75306
## 15 2023-03-01  96.44852  99.35893  96.86344
## 16 2023-04-01  95.57216  97.09987  95.76631
## 17 2023-05-01  96.07756  97.20263  95.25199
## 18 2023-06-01  96.01987  98.16253  96.68639
## 19 2023-07-01  98.06441 100.04565  97.70474
## 20 2023-08-01  98.71224 100.28066  98.70353
## 21 2023-09-01  97.60887  99.94061  97.89990
## 22 2023-10-01  97.69950  99.66011  98.34201
## 23 2023-11-01  99.70443 101.96982  99.87406
## 24 2023-12-01 100.64744 103.25350 101.62973

Ostateczne wyniki otrzymujemy biorąc średnią dla każdego miesiąca ze wszystkich metod:

3. Przewidywanie

Wyznaczona przez nas predykcja, łącząca w sobie modele ARIMA, Prophet oraz TBATS jako średnia arytmetyczna wszystkich prognoz prezentuje się następująco:

##       Date  Wskaznik
## 1  2022-01 103.52286
## 2  2022-02  99.45822
## 3  2022-03  97.03237
## 4  2022-04  95.89917
## 5  2022-05  95.93708
## 6  2022-06  96.67318
## 7  2022-07  98.27339
## 8  2022-08  98.94585
## 9  2022-09  98.25720
## 10 2022-10  98.38541
## 11 2022-11 100.30262
## 12 2022-12 101.61810
## 13 2023-01 104.21917
## 14 2023-02  99.97631
## 15 2023-03  97.55696
## 16 2023-04  96.14611
## 17 2023-05  96.17739
## 18 2023-06  96.95626
## 19 2023-07  98.60494
## 20 2023-08  99.23214
## 21 2023-09  98.48312
## 22 2023-10  98.56721
## 23 2023-11 100.51611
## 24 2023-12 101.84355

4. Podsumowanie

Literatura

  1. Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. J. (1990). STL: A seasonal-trend decomposition procedure based on loess. Journal of Official Statistics, 6(1), 3–33.

  2. regresja

  3. arima

  4. Taylor SJ, Letham B. Forecasting at Scale. The American Statistician. 2018

  5. De Livera A.M., Hyndman R.J., Snyder R.D., Forecasting time series with complex seasonal patterns using exponential smoothing. Journal of the American Statistical Association 2011, 106, 1513–1527.